#ifndef _SIUT_SIMD_STANDARD_BSPHERE_HPP_
#define _SIUT_SIMD_STANDARD_BSPHERE_HPP_

#include "Vec3f.hpp"
#include "BBox3f.hpp"

#ifdef DEBUG_BS
#define THROW_BS
#include <strstream>
#include <stdexcept>
#endif

#ifndef FLT_EPSILON
#define FLT_EPSILON 0.001
#endif

namespace siut
{
namespace simd
{

  struct BSphere3f
  {
    bool empty_;
    Vec3f origin_;
    float radius_;
    float radius_sqr_;

    BSphere3f()
      :empty_(true)
    {
      
    }

#ifndef NVCC
    BSphere3f(Vec3f orig, float rad)
      : empty_(false),
	origin_(orig),
	radius_(rad)	
    {
      radius_sqr_ = rad * rad;
#ifdef THROW_BS
      sanityCheck(this);
#endif
    }

    BSphere3f(const BBox3f &b)
      : empty_(false)
    {
      origin_ = lerp(getMin(b), getMax(b), 0.5f);
      radius_ = 0.5f * euclideanDistance(getMin(b), getMax(b));
      radius_sqr_ = radius_ * radius_;
    }    
#endif
  };  
  void clear(BSphere3f &bs);
  void updateDerived(BSphere3f &bs);
  bool containedClosed(const BSphere3f &bs, const Vec3f &p);
  const Vec3f& getOrigin(const BSphere3f &bs);
  float getRadius(const BSphere3f &bs);
  void include(BSphere3f &bs, const Vec3f &p);
  bool isEmpty(const BSphere3f &bs);
  void updateDerived(BSphere3f &bs);
#ifdef THROW_BS
  void sanityCheck(BSphere3f &bs);
#endif

  inline void clear(BSphere3f &bs)
  {
    bs.empty_ = true;
  }

  inline bool containedClosed(const BSphere3f &bs, const Vec3f &p)
  {
    if(bs.empty_)
      {
	return false;
      }

    return euclideanDistanceSquared(bs.origin_, p) <= bs.radius_sqr_;
  }

  inline const Vec3f& getOrigin(const BSphere3f &bs)
  {
    return bs.origin_;
  }

  inline float getRadius(const BSphere3f &bs)
  {
    return bs.radius_;
  }

  inline void include(BSphere3f &bs, const Vec3f &p)
  {
    if(bs.empty_)
      {
	bs.origin_ = p;
	bs.radius_ = FLT_EPSILON;
	bs.empty_ = false;
	updateDerived(bs);
      }
    else
      {
	float dSqrd = euclideanDistanceSquared(bs.origin_, p);
	if(dSqrd < bs.radius_sqr_)
	  {
	    return;
	  }
	bs.radius_ = std::sqrt(dSqrd);
	bs.radius_sqr_ = dSqrd;
      }
  }

  inline bool isEmpty(const BSphere3f &bs)
  {
    return bs.empty_;
  }

  inline void updateDerived(BSphere3f &bs)
  {
    bs.radius_sqr_ = bs.radius_ * bs.radius_;
  }
  
#ifdef THROW_BS
  void sanityCheck(BSphere3f &bs)
  {
    if(bs.empty_)
      {
	return;
      }
    if(bs.radius_ < 0.f)
      {
	stringstream s;
	s << "negative radius of bounding sphere in " << __FILE__ << " at " << __LINE__ << std::endl;
	throw runtime_error(s.str());
      }
  }
#endif
} //end namespaces
}


#endif
